library(Hmisc)
library(plyr)
library(latticeExtra)
library(ggplot2)
library(reshape2)
library(MASS)
library(lme4)
library(lmerTest)
library(simpleboot)
library(gridExtra)
library(cowplot)
library(coda)
library(runjags)

reorganizedGroups <- read.csv("RR2.csv")

#Fig 1
xyplot(distanceReference~round,groups=TREAT, data=reorganizedGroups, type="a",scales="free"
       ,ylab="Distance from the reference",xlab="Round",grid=T
       ,ylim = c(0.6,9)#, xlim = c(1.5,9.4)
       ,auto.key = list(text=c("RR0","RR1","RR2","RR3","RR4"), points=F, lines=T,corner = c(0, 0.9) )
       ,par.settings = (theme = custom.theme(symbol = brewer.pal(12, "Paired"), lwd = 3, lty = 1:5)))

## Bayesian parameter estimation of our model#bayesian estimation works best if controls are demeaned

# the spline against which we compare all other treatments
loessRr0 <- loess(distanceReference ~ round,data=subset(reorganizedGroups,reorganizedGroups$TREAT==1)) 


demean <- function(names,data)
{
  sapply(names,function(n) (data[[n]]<<-(data[[n]]-mean(data[[n]]))/sd(data[[n]])))
  data
}
#demeaning all variables were we are not interested in the cofficients
reorganizedGroups<-demean(names(reorganizedGroups[,c(7,8)]),reorganizedGroups)

initJags<-list()
initJags[[1]]<-list(.RNG.seed=1,.RNG.name="base::Mersenne-Twister")
initJags[[2]]<-list(.RNG.seed=2,.RNG.name="base::Super-Duper")
initJags[[3]]<-list(.RNG.seed=3,.RNG.name="base::Wichmann-Hill")
initJags[[4]]<-list(.RNG.seed=4,.RNG.name="lecuyer::RngStream") 
NumberOfChains=4;

##############
#For ESA
reorganizedGroups$rr <- (reorganizedGroups$TREAT - 1)/10
reorganizedGroups$rrSq <- reorganizedGroups$rr * reorganizedGroups$rr

#### ESA model with average spline

loessRr <- loess(distanceReference ~ round,data=reorganizedGroups)
reorganizedGroups$predictionAll <-predict(loessRr0,newdata = reorganizedGroups$round)

set.seed(123)
X.model <- 'model {
for (i in 1:length(distanceReference)) {
distanceReference[i] ~ dnorm(predictionAll[i] + betaRr * rr[i] + betaRrSq * rrSq[i] + betaPol * pol[i] + betaEnv * env[i] +
reGroup[Group[i]], tau1)
}
# random effects:
for (j in 1:max(Group)) {reGroup[j]   ~ dnorm(0,tau2)}

# fixed effects
betaRr ~ dnorm (0,.00001)
betaRrSq ~ dnorm (0,.00001)
betaPol ~ dnorm (0,.00001)
betaEnv ~ dnorm (0,.00001)
tau1     ~ dgamma(.01,.01)
tau2     ~ dgamma(.01,.01)
}'

X.jags<-run.jags(model=X.model,data=reorganizedGroups,burnin = 5000, sample = 10000,
                 monitor=c("betaRr","betaRrSq","betaPol","betaEnv"),
                 inits=initJags,method="parallel",n.chains=4)
X.df<-as.data.frame(as.mcmc(X.jags))
summaryJags <- summary(X.jags)

#gelman.diag(X.jags)
bayesX.mean <- mean(with(X.df,0<betaRr))
evidenceRr <-bayesX.mean
bayesFactorRr <- evidenceRr/(1-evidenceRr)

bayesX.mean <- mean(with(X.df,0>=betaRrSq))
evidenceRrSq <-bayesX.mean
bayesFactorRrSq <- evidenceRrSq/(1-evidenceRrSq)

bayesX.mean <- mean(with(X.df,betaPol<0))
evidencePol <-bayesX.mean
bayesFactorPol <- evidencePol/(1-evidencePol)

bayesX.mean <- mean(with(X.df,betaEnv>0))
evidenceEnv <-bayesX.mean
bayesFactorEnv <- evidenceEnv/(1-evidenceEnv)